Migration of a moonlet in a ring of solid particles : 
Theory and application to Saturn's propellers. 
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ABSTRACT 

Hundred meter sized objects have been identified by the Cassini spacecraft in Saturn's A ring 
through the so-caUed "propeller" features they create in the ring. These moonlets should migrate, 
due to their gravitational interaction with the ring; in fact, some orbital variation have been 
detected. The standard theory of type I migration of planets in protoplanetary disks can't be 
applied to the ring system, as it is pressureless. Thus, we compute the differential torque felt by 
a moonlet embedded in a two-dimensional disk of solid particles, with fiat surface density profile, 
both analytically and numerically. We find that the corresponding migration rate is too small to 
explain the observed variations of the propeller's orbit in Saturn's A-ring. 

However, local density fluctuations (due to gravity wakes in the marginally gravitationally 
stable A-ring) may exert a stochastic torque on a moonlet. Our simulations show that this torque 
can be large enough to account for the observations, depending on the parameters of the rings. 
We find that on time scales of several years the migration of propellers is likely to be dominated 
by stochastic effects (while the former, non-stochastic migration dominates after ~ 10^~^ years). 
In that case, the migration rates provided by observations so far suggests that the surface density 
of the A ring should be of the order of 700 kgm~^. The age of the propellers shouldn't exceed 1 
to 100 million years, depending on the dominant migration regime. 

Subject headings: planets and satellites: dynamical evolution and stability — planets and satellites: 
individual (Saturn) — planets and satellites: rings -L planet-disk interactions 



1. Introduction 

The theory of disk-satelHte interactions has 
for the most part been developed after Voy- 
ager's encounter with Saturn. The satelhtes 
that orbit beyond the outer edge of the rings 
perturb the dynamics of the particles compos- 
ing the rings. This leads to an exchange of 
angular momentum between the rings and the 
satellites, and t o the formation of d e nsity waves 
in the rings. Lin and Papaloizou ( 1979t ) and 



tudinal extent about 3 km ( Tiscareno et al. 20061 



Goldreich and Tremaind (|l979l Il980l ) have cal- 
culated, using two different methods, the total 
torque exerted by a satellite on a disk interior 
(or exterior) to its orbit. This torque is called 
the one-sided Lindblad torque, because it can be 
computed as the sum of the torques exerted at 
Lindblad resonances with the secondary body. In 
the lowest order, local approximation, the inner 
and outer torques are equal and opposite : the 
torque exerted by a satellite on a disk located 
inside its orbit is negative, with the same abso- 
lute value as the positive torque exerted on a disk 
located outside the orbit. 

Reciprocally, a disk exerts a torque on the sec- 
ondary body. When the strictly local approxi- 
mation is relaxed, the inner and outer torque s 
are not exactly equal and opposite ( Wardl ll986'l. 
Their sum, called the differential Lindblad torque, 
is generally negative. As a consequence, the or- 
bital angular momentum of a body embedded in 
a disk decreases, and so does its semi-major axis 
(on circular Keplerian orbits, the orbital angular 
momentum is proportional to the square root of 
the semi- major axis). This is planetary migration 
of type I (|Wardlll997l) . So far, this phenomenon 
has been mainly studied in the frame of plan- 
ets embedded in protop lanetary gaseous disks (see 
Papaloizou et al.ll2007l for a review). 



The Cassini spacecraft has been orbiting the 
Saturnian ring system since 2004, offering the 
possibility to observe the coupled evolution of the 
ring system and the satellites. Due to short orbital 
timescales (1 year is equivalent to about 700 orbits 
of the A ring) it may be possible to observe the ex- 
change of angular momentum between the two sys- 
tems. One of the most striking discoveries of the 
Cassini spacecraft is the observation of propeller 
shaped features in the A ring (located between 
122 000 and 137 000 km from Saturn), with longi- 



I2OO8I : ISremcevic et allbOOTt ). They are most prob 
ably caused by the presence of moonlets about 
hundred meters in size, emb edded in the ring, and 
scattering ring particles (jSpahn and Sremcevic 
:2000i ). As they are embedded in the ring, these 
small bodies should exchange angular momentum 



with the ring, and migrate ( Crida et al. 20091 ) 



This migration could be detected by Cassini ob- 
servations through the cumulative lag, or advance 
with time t of the orbital longitude 4> induced by 
a small variation of the semi-major axis and the 
angular velocity ft {54) = 5^ y. 6t), offering for 
the first time the possibility to confront directly 
the planetary migration theory with observations, 
and to give insights and constrains on the physical 
properties of the rings and of the moonlets. 

In this paper, we address the question of the 
theoretical migration rate of these propellers, us- 
ing both numerical and analytical approaches. 
The theory is then confronted to observations. In 
Sect. [21 we review the standard theory of type I 
migration ; the differences between migration in 
protoplanetary disks and in Saturn's rings are ex- 
plained, showing the need for a new calculation of 
the migration rate of embedded moonlets. This 
rate is given in Sect. [3] for an homogeneous, ax- 
isymmetric disk with a flat surface density profile, 
as a result of numerical computation in Sect. 13.11 
and analytical calculation in Sect. 13.21 In Sect.HJ 
we consider the effect of density fluctuations in the 
rings, in particular the role of short-lived gravitat- 
ing clumps, also called gravity wa kes, which are 



known to be numerous in the A ring (jColwell et al 



20061 ). We then conclude in Sect. [5] on what our 
model tells us on the properties of the rings, given 
the observed migration rates. 

2. Review of Type I migration and the dif- 
ferential Lindblad torque 

In protoplanetary gaseous disks, the perturba- 
tion caused by a terrestrial planet leads to the for- 
mation of a one armed spiral density wave, leading 
the planet in the inner disk, and trailing behind 
the planet in the outer disk. This wave is pressure 
supported and generally called the wake, but it has 
nothing to do with the gravity wakes mentioned 
above : the latter are local features, while the 
planet wake spirals through the whole disk. The 
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planet wake carries angular momentum, so that 
the angular momentum given b y the planet to the 
disk is not deposited locally (e.g. lCrida et al.lliood . 
Appendix C). Therefore, the disk profile is hardly 
modified in this linear regime. Still, the negative 
torque exerted by the outer disk on the planet 
through the wake is larger in absolute value than 
the positive torque from the inner disk. Without 
going int o the details (for which the reader is rcf- 
ereed to I WardI HqOtI ) . the main reason the outer 
disk wins over the inner disk lies in pressure ef- 
fects : from the dispersion equation of a pressure 
supported wave, one finds that the location VL^rn 
where the wave with azimuthal mode number m, 
corresponding to the mth Lindblad resonance with 
the planet, is launched, is not exactly the location 
of the resonance given by Kepler's laws. The shift 
is not symmetrical with respect to the planet posi- 
tion for inner and outer resonances, but favors the 
outer ones. As a consequence, the planet feels a 
negative total torque, ca lled the differential L ind- 
blad toque, and given bv lTanaka et all (|2002l ) : 



Tdiff = -Cq Y,rp Qp h 



2,-2 



(1) 



where the index p refers to the planet, rp being 
the radius of its orbit and flp its angular velocity, 
q is the planet to primary mass ratio, S is the 
surface density of the disk in the neighborhood 
of the planetary orbit. Finally, h is the aspect 
ratio of the disk, being the ratio between its scale 
height and the distance to the central body, being 
proportional to the square root of the gas pressure. 
In a typical protoplanetary disk, h 0.05. The 
numerical coefficient C is given by C = 2.340 — 
0.099i^, where ^ is the index of the power law of 
the density profile: E cx r~^. 

This result is robust. In particular, the value of 
the negative torque is almost independent on the 
slope of the density profile This is due to the so- 
called pressure buffer: the resona nces a.re shi fted 
when the density gradient varies (|Wardlll997() . If 



the disk were pressureless, then the expression of 
C would be completely different. Also, the aspect 
ratio h in Eq. ^ appears because r^.m doesn't 
converge towards rp when m tends to infinity, but 
towards rp(l ± 2/i/3), due to pressure effects. To 
sum up, the gas pressure plays a fundamental role 
in type I migration. 

It should be mentioned for completeness that, 
in addition to the differential Lindblad torque 



discussed above, the horseshoe drag - exerted on 
the planet by the gas on horseshoe orbits around 
the planetary orbit - pl ays a signi ficant role in 
type I migration (see e.g. WardlfigQl : MassefeoOl 



iBaruteau and Masset"2008'; 'Kley and Crida"200g ■ 



Paardekoopcr a nd Papaloizou.2009;.Paardekooper et al, 
fe009h . 

In contrast to gaseous protoplanetary disks, 
pressure effects in Saturn's rings are not impor- 
tant. The aspect ratio h is of the order of 10^^. 
The spiral density waves that are observed in the 
A ring are gravity supported, not pressure sup- 
ported. Thus, the standard theory of type I migra- 
tion does not apply. In particular, the spiral planet 
wake doesn't appear. The interaction of the moon- 
let responsible for the propeller structure with the 
disk is observed to take place within a few hundred 
kilometers. Resonances with m > 10'^ are located 
within this distance and should play a significant 
role. However, in the standard type I migration 
the important resonances have m ~ ^ lO*". 
Therefore, Eq. ([!]) can't be directly applied to a 
moonlet in Saturn's rings. A new approach is 
needed, adapted to the two main characteristics 
of the problem : the fact that the rings are made 
of solid particles, and the fact that the interaction 
is taking place very close to the moonlet. 

3. The ring-moonlet interaction 

In this section, we compute the interaction be- 
tween a moonlet and a ring test particle. In this 
analysis, the gravity of the other ring particles is 
neglected. This leads to the torque exerted on 
the moonlet by an initially unperturbed, homoge- 
neous ring. In subsection 13.11 the computation is 
performed numerically. In subsection [321 it is de- 
rived analytically. Both results are in agreement, 
and a corresponding migration rate for the moon- 
let is given and discussed in subsection 13.31 

From now on, m denotes the mass of the moon- 
let (and not anymore the order of a resonance). 
The moonlet, has a circular orbit of radius r„j 
around the central planet of mass M. The grav- 
itational potential due to the moonlet is The 
radial and azimuthal components of the equation 
of motion of a ring particle in two dimensional 
cylindrical polar coordinates (r, 0) are 



dt^ 



d4) 



a* GM 

dr r^ 



(2) 
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and 



(3) 



dt^ '^^[dtj [jdt J ~~r~d^' 

The angular velocity of the moonlet is a; = 
sjGM jrrr? ■ Let rg be the radius of the initially 
circular orbit of a test particle and i7 = \fGMjr^ 
its angular velocity. We note & = — is the 
impact parameter, and h — b/rn the normalized 



impact parameter, where th = Vm {jm) 
Hill radius of the moonlet. 



1/3 . 



is the 



3.1. Numerical computation of the ring 
moonlet interaction 

In this section, the numerical integration of the 
above equations of motion is performed, in order 
to find the trajectories of ring particles in the pres- 
ence of a perturbing moonlet in the frame corotat- 
ing with the moonlet. To measure the tiny asym- 
metry between the inner and the outer part of the 
ring, the full equations are integrated, without lin- 
earizatioi i or simplification . A Bulirsch-Stoer al- 
gorithm ( Press et al. 19921 ) is used, and a Taylor 
expansion is performed in the code when neces- 
sary to subtract accurately large numbers, in order 
to achieve machine double precision (10^^^). We 
have checked that the Jacobi constant is conserved 
to this precision along the trajectories. Examples 
of obtained trajectories are given in Fig. [TJ 

It is well known within the framework of the re- 
stricted 3-body problem that if |6| is small enough, 
the test particle has a horseshoe shaped orbit, 
while if |6| is larger than ^ 2.5, the test particle is 
circulating, and scattered into an eccentric orbit. 
This can be seen in Fig. [T] We perform many nu- 
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Fig. 1. — Trajectories of test particles perturbed 
by a moonlet of mass m — 3 x 10^^^ M, located 
at (r = rm,4> = 0) (that is at (0,1) in the plot), 
in the frame corotating with the moonlet. Dashed 
circle : orbit of the moonlet. 



merical integrations with various b in the case of 
a moonlet of mass to = 3 x 10~^'^M, starting the 
particle at an azimuth \(f)Q\ = 3000 rj^/r^ = 0.3 
(where the moonlet is at = 0). This angle is 
large enough so that at this location the influence 
of the moonlet is negligible and the orbital pa- 
rameters of the test particle are not disturbed, as 
will be checked later. We find that the horseshoe 
regime occurs for b < 1.8, and the scattered regime 
occurs for b > 2.5. For 1.774 < 6 < 2.503, how- 
ever, the trajectory approaches the center of the 
moonlet to within a distance smaller than 0.95 rjy. 
In that case, if one assumes the moonlet is a point 
mass, the test particle eventually leaves the Hill 
sphere, either on a horseshoe or a circulating tra- 
jectory, but the outcome changes several times 
with increasing b. In the case we are concerned 
about here, the moonlet most likely almost fills 
its Roche lobe, and therefore we stop the inte- 
gration of the trajectory as soon as the distance 
between the test particle and the moonlet is less 
than 0.95 r//, assuming a collision. 

The specific orbital angular momentum J = 
r'^{dip/dt) of the test particles is computed along 
the trajectories. Angular momentum is exchanged 
during the close encounter with the moonlet. For 
b ^ 2.503, the test particle is scattered onto an 
eccentric orbit of larger angular momentum than 
initially, which results in a gain in angular mo- 
mentum. The variation of orbital angular momen- 
tum along the trajectory is shown in the bottom 
panel of Fig. [2] for the case b — 3, where the top 
panel is the trajectory. The difference in angu- 
lar momentum between the initial circular orbit 
at (j}o = 0.3sgn(b) and the end of the integration, 
when = 0.3 again, is noted A J. In the figure, 
only the interval —0.05 < < 0.05 is displayed, 
for convenience. Most of the exchange of angular 
momentum occurs when < 0.01. 

Figure [3] shows | A J| (top thick curve) as a func- 
tion of b, in units of the specific angular momen- 
tum of the moonlet Jm = f'^^- For < 6 < 1.774, 
the horseshoe trajectory corresponds to a U-turn 
towards the central planet, and to a loss of angular 
momentum for the test particle. More precisely, as 
for circular orbits J cx r^/^, one expects for such a 
U-turn A J/ J = i Ar — — b/r^ ; this is indeed the 
case for b < 1.3. In the case where the test particle 
collides with the moonlet, we assume that it gives 
all its orbital angular momentum to the moon- 
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let: AJ = rg^n - J,„, so that AJ/J sa b/{2r^). 
This also appears in Fig. [3l The opposite holds 
for 6 < 0. 

Computing AJ as a function of 5 to numer- 
ical precision enables us to also compute the 
difference between the inner and outer disk : 
SJ{b) = AJ{b) + AJ(-6). This quantity is 
small with respect to AJ{b), but nonetheless well 
determined and converged in our simulations : 
AJ(6) + AJ(— 6) is constant after the encounter 
to a precision better than 0.5% for all \<j)\ > 0.02. 
This validates our choice of (pQ. In Fig. [3l the 
bottom thick dashed curve shows dJ in the same 
scale as |AJ|. We see that SJ > for all 6 > 
and that SJ <C A J, with 

(5 J/A J « 5 X lO"-* 6 (4) 

for circulating trajectories, and 

6J/\AJ\ w 1.17 X IQ-'^b 

for horseshoe orbits. In the following subsection, 
the empirically found Eq. Q is derived analyti- 
cally and justified. 

3.2. Analytic model for the ring moonlet 
interaction 

In this section, we consider only circulating tra- 
jectories. Developing to second order the exchange 
of angular momentum during an encounter with 
the moonlet A J, we can find the asymmetry SJ. 

3.2.1. Solution for the perturbed moonlet orbit 

Let us start again from Eq. ^ and ([3]). The 
ring particle is assumed to be on an unperturbed 
circular orbit of radius tq = + b. It orbits with 
angular velocity Vt = GM/r^^ such that = ili, 
where without loss of generality we have defined 
the origin of time t — Q to he when the particle 
is at = 0. Under the perturbation induced by 
^, the particle moves to r — r^ + x, and cj) ~ 
Clt+y/Vm., where x and y are assumed to be small. 
Linearizing Eq. ^ and ([3]) about the circular orbit 
state, we obtain equations for x and y in the form 



— - 2n-f- - sn^x 

dt^ dt 



dt^ ^ dt 



dr 

1 a* 

ro dcj) 



and 



(5) 
(6) 
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Fig. 2. — Top panel : trajectory of a test particle 
with impact parameter b = Z \ the motion of the 
particle is toward negative <j}. Bottom panel : vari- 
ation of the specific orbital angular momentum J 
of the same particle along its trajectory. 
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Fig. 3. — Angular momentum exchanges during 
one close encounter, as a function of the impact pa- 
rameter. Top, thick, red curve: AJ(6), from nu- 
merical simulations. Green, thin, dashed, straight 
line: AJ(&), as given by Eq. ([5^ . Bottom, thick, 
dark blue, long-dashed curve: SJ{b), from nu- 
merical simulations. Thin, light blue, dash-dotted 
line: 5J{b) as given by Eq. ([M)) . taking A J from 
the simulations. Orange, thin, double- and triple- 
dashed lines: (b/rm)Jm, and {b/2rm)Jm, respec- 
tively, to compare with |AJ|. 
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Here the subscript denotes evaluation on the un- 
perturbed particle orbit. 

We suppose that the perturbation is induced by 
a moonlet of mass m that is on a circular orbit of 
radius and has an angular velocity uj. Then its 
azimuthal coordinate 0™ = ujt + 4>rn,Q, with 0m, o 
being a constant. The perturbing potential 

Gm 



V To + '"m ~ 2r,„ro cos((/i ~ 



(7) 



becomes a function of time through substituting 
— 4>m = {^l — uj)t — (jimfi therein. 
Thus we have 



1 

ro d(j) 



1 



ro(il — ll)) dt 



(8) 



Using this in Eq. (|6|) and integrating with respect 
to time, we obtain 



1 + 2^^-- <n ^ 

at ro[il — uj) 



(9) 



which when combined with Eq. ([5]) gives an equa- 
tion for X in the form 



fx 



n^x 



dr 



2m 



ro(rj - uj) 



S (10) 



3.2.2. Solution of the linearized equations 

To solve Eq. PU)) . we note that the perturbing 
potential Eq. ([7]) evaluated on the unperturbed 
orbits is a periodic function of time with period 
27r/|a;— r2| = 2tt//3. Thus we should look for a peri- 
odic response. In order to do this we have to intro- 
duce a small frictional term into Eq. ([TU]) to enable 
transients to decay and a net torque on the moon- 
let to be set up. When the frictional term is small 
it is expected that the res ulting torq ue should not 
depend on it (e.g. Goldre ich and Tr cmainc 1980). 
Hence we add a frictional term ^(dx/dt) to the 
left hand side of Eq. ([TU]), where j/il is a small 
constant parameter so that it now reads 



d^x dx 
'd^^^'dt 



9* 

dr 



ro(fi - uj) 



(11) 

As the potential is periodic in time we can adopt 
a Fourier series of the form 



S = 



n— — oo 



(12) 



where it is implicit that the real parts of such com- 
plex expressions is to be taken, and 

Sn = T- S{t) exp{-inl3t)dt . (13) 

27r Jo 

The periodic solution of Eq. (1111) is now readily 
written down as 



E 



n— — oo 



Sn exp(m/3t) 
(rj2 - „2^2 _^ ^^nl3) ■ 



(14) 



We may write this in terms of a Green's function 
defined through 



G(r) 



1 

2^ 



n— OO 

E 

n— — oo 



exp(in/3r) 



Then the solution for x may be written 



x^ P 



27r//3 



S{t-t')G{t')dt'. 



(15) 



(16) 



Note that as the orbit of a ring particle relative 
to the planet is periodic, the solution given by 
Eq. (fTB)) includes the effects of infinite numbers of 
repeating encounters. However, we wish to con- 
sider the case when dissipative effects, although 
weak, are strong enough to recircularize orbits be- 
tween encounters in which case they will be in- 
dependent of each other. This condition requires 
that 7/|a; — f7| = 7//? ^ 1. This is equivalent to re- 
quiring that the damping time scale be short com- 
pared to the relative orbital period between moon- 
let and ring particle. On account of the length 
scale of the encounters of interest being compara- 
ble to the Hill radius of the moonlet, this is much 
longer than the orbital period itself, so that we 
may adopt the ordering 



7/|a;-f7|=7//3>l>7/17. 



(17) 



In order to make use of the above ordering we write 
down the form of the Green's function derived in 
the appendix (see Eq. (jA7l) ) valid for < t < 
2^//3. 

G(r) = 



sin(a;^T) 



-T'^/'^ sin(a;^(r-27r/^)) 



uj-iP [I + e-27'r//3 _ 2e-T^/'3 cos{2tt Lo ^ / [3)] 



(18) 
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where = -y/fi^ _ ^2^4_ ^j^g function is de- 
fined elsewhere through its periodicity with pe- 
riod 27r//3. Making use of the inequahty Eq. (fTTj) 
we may replace the Green's function Eq. by 
the simple expression 



G(r) 



exp(— 7r/2) sin(a;^T) 



(19) 



Then the solution Eq. gives 



X = A(i)e"'^*/2 sin(w^t) + B{t)e-^^''^ cos{u:^t), 

(20) 

where 



= — / 5'(t') exp(7i72) cos{uj^t')dt' 



(21) 



and 



B{t) = / S{t') exp(7<72) sin{uj-ft')dt' . 

(22) 

To make use of the above expressions, we con- 
sider the situation when the ring particle has a 
close encounter with the moonlet at time t — 0, 
thus we take (j)m.o ~ (we note that a non zero 
'/'m.o can be dealt with by rotating the coordi- 
nate system and shifting the origin of time). The 
source term S is then expected to be highly peaked 
around t' = 0, and almost all of the contributions 
to the above integrals will occur for \t'\ <~ 27r/57. 
Furthermore, during this dynamical interaction, 
dissipation will be negligible. Thus, if we are in- 
terested in times after the main interaction, but 
before significant dissipation takes place, we may 
set 7 = and extend the limits of the integra- 
tion to ±cxD. However, in practice one may have to 
apply a cut off to the potential at large distances 
from the moonlet in order to do that (see below). 
But this should not matter if the important inter- 
action occurs when the moonlet and ring particle 
are close. 

Then we simply have 



A = 



s{t') cos{nt')dt' 



(23) 



and 



S{t')sm{nt')dt'. (24) 



Thus A and B are constants representing epicyclic 
oscillation amplitudes induced after the close ap- 
proach of the ring particle to the moonlet. 

We remark that the approximations made in 
obtaining Eq. ((23)) and Eq. (l24t relate to how dis- 
sipation is treated. There has been no assump- 
tion that the particle trajectories are symmetric 
on opposite sides of the moonlet so that curvature 
effects remain fully incorporated during particle 
moonlet encounters. When dissipation is negligi- 
ble during the encounter, then immediately after- 
ward an epicyclic oscillation is established. The 
assumption that dissipation circularizes orbits be- 
tween encounters implies that we should consider 
approaching ring particles to be on circular orbits. 
The above discussion indicates that errors associ- 
ated with this assumption are exponentially small. 

3.2.3. Angular momentum transfer 

For the set up considered here, symmetry con- 
siderations imply that S{t) is an even function of 
time (see below), so that B = Q. The generation 
of the epicyclic oscillation is associated with an 
angular momentum transfer between the moonlet 
and particle. To find this we firstly note that 



A J = VGM i y^a/(l -e2) 



(25) 



where a / and e are the post encounter semi- major 
axis and eccentricity of the particle. We also note 
that the Jacobi constant implies that the change of 
the particle orbital energy and angular momentum 
are related by l\E = GM[l/(2ro) - l/(2a/)] = 
Lo A J. This can be used to eliminate a/ in Eq. (P5|) 
after which AJ may be found correct to second 
order in e = A/r^ with the result that AJ = 
ri^^^/[2(a; — O)]. This in turn may be simply de- 
termined after evaluating A. Note that AJ < 
for particles interior to the moonlet which have 
Q. > oj and conversely A J > for particles orbit- 
ing exterior to the moonlet. 

3.2.4. Development of the perturbing potential 
We now consider 



dr 



ro(f^-w) 



(26) 
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We begin by recalling that 



Gm 



Gm 



. ■ (27) 

V ^0 + - 2ror,m cos(/3t) 

In order to evaluate the Fourier transform as spec- 
ified by Eq. ([23|). which was derived under the 
assumption that the interaction occurs only near 
closest approach, we must truncate the potential 
at large \t\. As the encounter takes place over a 
time <C this can be achieved by replacing 

cos(/3i) in Eq. dS?]) by 1 - ii'^t^ 12. Note that a 
dimensionless estimate of the error involved is of 
order (JijiS)^ ~ {th /rmf' ^ where rn is the Hill ra- 
dius of the moonlet. This is small enough that the 
leading order asymmetry in the angular momen- 
tum transferred to orbits with the same impact 
parameter on either side of the disk can be esti- 
mated. 

As the first stage in evaluating the Fourier 
transform of S specified in Eq. (|23|) that gives the 
epicyclic amplitude, we evaluate 

1 

C ^ - *cos(f7t)di 



Gm cos(0) 



This can also be expressed as 

2GmKo{^o) 



C 



(28) 



(29) 



where = {^Vq - rm\) / {P\/nyr^), and K-j de- 
notes the modified Bessel function of the second 
kind of order j. 

3.2.5. Total angular momentum exchange 

We may now use the above expression together 
with Eq. (pS)) to evaluate the epicyclic amplitude 
Eq. ([25]) (noting that the radial derivative is with 
respect to rg with other quantities held fixed) so 
obtaining 



A 



2Gm 



1 



2n 



^i(eo)eo 



(30) 



The associated angular momentum exchanged is 
then given by 



2(Gm)2 



' Q' m 



{UJ - f7)3 



2n 



(O-w) 



-ifi(eo)Co 



(r-0 



(31) 



In a strictly local approximation under which 
the inner and outer sides are symmetric, the con- 
tributions from orbits equidistant from the moon- 
let would cancel, leaving the net result to be de- 
termined by the surface density profile. However, 
although we have assumed the interactions are lo- 
cal, we did not assume symmetry between the ex- 
terior and interior orbits. Accordingly we evalu- 
ate the difference in the magnitude of AJ eval- 
uated from orbits equidistant from the moonlet : 
''o — Tmib. The leading order contribution to A J 
is symmetric in b. The lowest order contribution 
is antisymmetric and accordingly leads to cancel- 
lation between the two sides. We make use of the 
expansions = 2/3 — b/{2rm) + 0{{b/rm)'^), and 
2n/in - w) = -4r„/(36)(l - 6/(4r„)) + 0{b/r^) 
together with standard properties of Bessel func- 
tions to write 



(32) 



where 



3 (6i.,(2/3) + 3i.o(2/3))^ (33) 

4 (4ifo(2/3) + 2i^i(2/3) ^ ' 

Th e first order term of Eq. (1321) wa s already given 
by iGoldreich and Tremaine It is plotted 

as a straight green dashed line in Fig. [3] Our ex- 
pansion to second order enables us to go further, 
and to give the expression of the magnitude of the 
asymmetry between the two sides of the disk : 



SJ 
AJ 



2a\b\/r^ = 4.92|6|/r„ 



(34) 



It is such that for an orbit with a given impact 
parameter, the angular momentum exchanged in 
the outer disk is the larger. 

In the case studied numerically, we had rn = 
10-^, so that \b\/rm = 10^*b. Then, Eq. ^ 
remarkably agrees with the numerical fit Eq. 
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The light blue dot-dashed curve in Fig. [3] displays 
4.92 X IO-^SAJ. 

In the context of the above, we note that ap- 
proximations made in obtaining equation Eq. pip 
such as effectively starting and truncating the in- 
teraction at some finite though large distance from 
the moonlet could conceivably lead to changes 
comparable to those given by Eq. (|34|) . However, 
such changes are again approximately symmetric 
for trajectories on both sides of the moonlet and 
thus approximately cancel so we do not expect 
such effects to significantly alter Eq. p4|) . 



3.3. Migration rate and discussion 

If the surface density of ring particles is S, the 
total rate of angular momentum transferred to the 
moonlet is 



dJ_ 



E AJ- 



disk 



2tt 



-dr rc 



(35) 



where the integral is taken over the disk. The 
particles exterior to the moonlet contribute neg- 
atively while those interior contribute positively. 
The cumulative torque exerted by the moonlet on 
the region of the ring located within a distance b 
to its orbit reads then : 

T^{b)= I T.{rra + h'){/\J{b'))\uj~n\db' (36) 

J -b 

The normalized cumulative torque 

TMI [(m/M)4/3(E/A/r„f )^ 

is plotted in Fig. 21 The proportionality to S is 
obvious; that Tc oc (m/Af)'*/^ is numerically ver- 
ified for 3 X 10^15 sC m/Af 3 X I Q-^ and has 



been already found analytically by IWardI (Il99lh 
for the horseshoe drag in a similar context. 

Most of the total torque comes from scattered, 
circulating particles, in particular the ones with 
smallest impact parameter b sa 2.5. This makes 
the total torque sensitive to the physical size of the 
moonlet (taken as 0.95 r// here), as some particles 
colliding with the moonlet could be circulating if 
it were smaller. 

The role of the horseshoe drag appears to be 
non negligible, amounting to 
- 4.1( i:/Mr-^)(m/ Af)^/3 Afr,>^ The expres- 
sion of WardI ( 199l[ ) for the torque arising from 



material executing horseshoe turns, called the 
horseshoe drag, is for a Keplerian disk with flat 
density profile : 



(37) 



where w is the half-width of the horseshoe region. 
In our case, w = 1.774 rjj, which gives Ths = 
2.6(E/Afr-2)(m/A/")''/3 Mr^cu'^. The agreement 
is good because Ward's analysis is based only on 
geometrical effects and angular momentum vari- 
ation in a Keplerian disk, without any pressure 
effect. Therefore, it also applies in Saturn's ring. 
We remark that taking w = 2rH in Eq. (|37l) gives 
a perfect match with what we find numerically for 
the total horseshoe drag. 

In conclusion, from Fig. |4l the total torque felt 
by a moonlet of mass m on a circular orbit of ra- 
dius Tm around a planet of mass M can be written 
as 



T 



-17.8 



\Mrn 



(-) 



4/3 



Arr„V- (38) 



Note that to get the same dependency of the type I 
torque in the parameters of the system, one has 
to assume h cx TH/rm in Eq. ([ij ; however, in a 
protoplanetary disk, h is fixed and independent 
of the mass of the secondary body, so that this 
proportionality would not be justified. 

The torque is related to the migration speed 
through T — 0.5 mrm^(drm/dt). Hence we de- 
duce that 



dVyji 

~dt 



= -35.6- 



M 



'mNi/3 



(39) 



The migration rate is here proportional to the 
mass of the moonlet to the power 1/3, in con- 
trast to standard type I migration where drm/dt cx 
m. A numerical application to the case of an 
TO = 10-^^Afsaturn = 5-68 X 10® kg mooulct in 
orbit in the A-ring of density E = 400 kgm~^ 
at Tm = 130 000 km from Saturn gives drm/dt = 
—0.23 myr^^. Increasing the mass by two orders 
of magnitude to correspond to a radius of ^ 200 m 
speeds up the migration rate by a factor of only 
4.5 to ~ —1 niyr~^. 

After time t, a migrating propeller will be 
shifted longitudinally with respect to a corre- 
sponding non migrating one by a distance r„i Scjj — 
3ri |(ir,„/(ii| t^/4. For the above parameters, this 
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gives 713 year)]^ m. A shift of this magni- 
tude is potentiahy de tectable on a timescale of a 
year to a few years (IPorco et al ][2004 and note 
also Burns etall |20Q9| )^ Actually, migration of 



propellers has already been d etected () Burns et al 



20091: iTiscareno etaU l2010t ). During one time 



period of nearly a year, a particular propeller 
has been seen moving outward at a rate of ~ 
110 myr~^ ; and during a later similar time pe- 
riod, the same propeller has been seen moving in 
ward at a rate of ~ 40 myr~^ (ITiscareno et al 



201fll and personal communication) 

These observations are not compatible with the 
above theory. But we recall that the process of 
migration of a moonlet described above assumed a 
smooth particle disk with constant surface density. 
Here we note that there are features and mech- 
anisms that might produce a significantly faster 
migration rate, possibly in both directions inward 
and outward, and non constant in time. One can 
first think of a radial density gradient : as there is 
no pressure buffer here, this would directly affect 
the balance between the torques from the inner 
and outer parts of the ring. This would also affect 
the torque from the horseshoe region, which could 
turn positive. However, if the migration is gov- 
erned by the gradient of some quantity, it seems 
likely that the moonlet would have approached an 
extremum in that quantity, and thus should have 
attained a migration rate comparable to that esti- 
mated for a constant surface density. 

Another possibility resulting in the moonlet 
migrating faster than what the previous calcu- 
lation indicates, and possibly outward, is a run- 
away migration in a plan etesimals disk (jlda et al 



20001: iLevison et al.ll2007l for a review), similar to 
the type HI migration of planets in pr otoplane- 
tary disks ( Masset and Papaloizou I2OO3I ) . In this 
regime, the migration of the moonlet in the disk 
leads to a positive feedback on its migration rate, 
because of the material of the inner (resp. outer) 
disk making horseshoe U-turns to the outer (resp. 
inner) disk. This speeds up the migration, pos- 
sibly leading to a runaway. However, this leads 
inevitably to an asymmetry in the horseshoe re- 
gion, while the propeller structures observed are 
rather symmetrical. 

Finally, the A-ring of Saturn is not homoge- 
neous. It is close to gravitational instability, which 
should lead to the formation of gravity wakes and 



density fluctuations. The effect of these density 
fluctuations on the moonlet is studied in next sec- 
tion. 

4. The role of density fluctuations and re- 
sulting stochastic migration 

The analytic calculations and numerical simu- 
lations in the previous chapters assume an inflow 
of particles on circular orbits only perturbed by 
the nearby moonlet. However, we know that Sat- 
urn's A ring is mar ginally gravitationally stable 
('Daisaka et al. 2001). The Toomre Q parameter 
(,Toomre,196l . which is a measure of the impor- 
tance of self-gravity, is expected to be of the order 
of 2 ^ 7, indicating that the ring particles' mutual 
gravity is indeed a strong effect. It leads to the 
regular formation and dispersion of gravity wakes, 
which are local density enhancements elongated in 
parallel directions by the Keplerian shear. Those 
over-densities give rise to stochastic forces which 
act on the embedded moonlet. 

A very similar effect is expected to occur in 
protoplanetary disks. These disks are thought to 
be turbulent d ue to the magneto-rotatio nal in- 
stability (MRI, [Saibus^ndH^ie^ llMll) . The 
turbulent fluctuations create over-densities which 
interact gravitationally with embedded small 
mass planets. The stochastic forces make the 
planet undergo a random walk. An analytic 
model of this random wal k has been derived by 
iRein and Papaloizou (2009|). In the following, we 
apply this model to moonlets embedded in Sat- 
urn's rings. To do that, we need to get an estimate 
of the amplitude of the stochastic forces. 

4.1. Numerical calculations 



We perform three-dimensional simulatio ns of 
ring p articles, in a shearing box, similarly to Said 
(jl995l ). The simulations are done in a local cube 
of size H with shear periodic boundary conditions, 
and the origin of the box is fixed at a semi ma- 
i or axis of a = 130 O OP km. A BH tree code 
( Barnes and Hut 1986 ^ is used to calculate the 
self-gravity between ring particles and resolve in- 
elastic collisions. Collisions between particles are 
resolved using the instantaneous collision model 
and a velocity dependent coefficient of restitution 
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given bv lBridges et al. 



e{v) = min <^ 0.34 x 



( " ) 
\ 1 cm.s 1 / 



-0.234 



, 1 



(40) 



where v is the impact speed projected on the vec 
tor joining the centers of the two par ticles. The 
code is described in more detail in iRein et aT 
(|201fl[ ). 



The size of ring particles is not well constrained. 
Therefore and to be able to scale to different loca- 
tions in Saturn's rings, we perform multiple sim- 
ulations. For a given simulation, all the particles 
are assumed to be spherical and have the same size 
(or radius), which varies from simulation to simu- 
lation from 0.52 to 13 meters. The simulation pa- 
rameters are listed in Table [1] The nomenclature 
and physical pa rameters are, for ea s y corn parison, 
the same as in Lewis and Stewart ( 20091 ). as our 
simulations are similar to theirs. 

The moonlet is not taken into account in the 
simulations. We measure the specific gravitational 
force f (or acceleration) felt by a passive test- 
particle sitting at the origin. We calculate the 
force in two different ways, in order to avoid the 
singularity at the origin and to account for the 
physical size of the moonlet. In the first case, we 
use a cut off at the moonlet 's radius d and exclude 
all particles within that radius from the force cal- 
culation. In the second case we use a smoothed 
gravitational force per unit mass in the form 



f = 



Gm 



part 



d2 



(41) 



where f is the vector linking the origin to the 
particle and mpart is the mass of the particle. 
The smoothing length d is set equal to the moon- 
let's size. In a self-consistent simulation, one 
should include the moonlet with it's real physi- 
cal size. However, this goes beyond the scope of 
this paper and will be cons idered in future work 
( Rein and Papaloizoull20ld) . Our purpose here is 
to estimate the underlying stochastic fiuctuations 
in the migration rate that occur independently of 
the moonlet. This procedure is reasonable as long 
as the moonlet is in a steady state, namely if it 
doesn't accumulate or lose a large amount of mass 
over one orbit. In all our simulations we assume a 
moonlet size oi d = 200 m. 




Fig. 4. — Cumulative torque given by Eq. (|36|) . 
exerted by a moonlet on the region of the ring 
rm - b < r < r„i + h. 
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Fig. 5. — Azimuthal component of the specific 
gravitational force felt by the moonlet in ms"^ for 
simulations L2, S3 and SI (from top to bottom). 
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Name 




T 




Pp 




Tip 


H 




N 


L2 


13 m 


0.1 


0.5 


gcm" 


-'A 


885 kgm-^ 


5000 


m 


4808 


SI 


0.52 m 


0.1 


0.7 


gcm~ 


-3 


49.7 kgm"2 


1000 


m 


120 548 


S3 


1.3 m 


0.2 


0.7 


gem" 


-3 


246 kgm^2 


1000 


m 


38188 



Table 1: Simulation para meters. The first column identifies the simulation, following the convention of 
Lewis and Stewart (l2009l) . The second and third column give the size (or radius) of the particles and their 
density. The fourth and fifth column give the surface density and the size of the computational domain, 
respectively. The last column lists the number of particles. 



4.2. Results 

We measure the amplitude and the correlation 
time of the stochastic forces in all simulations. 
The results are listed in Table [H We also plot 
the time evolution of the azimuthal (y) force com- 
ponent in Fig. [5l Whereas the correlation time in 
all simulations is almost the same, the amplitude 
of the stochastic fluctuations varies by almost a 
factor of 10'^. The forces in the vertical direction 
are negligible and not presented here. The diffu- 
sion coefficient, being a measure of the strength 
of stochastic forces, is defin ed as _D = 2(/^)r (see 



Rein and Papaloizou 2009[ ). where (/^)^/^ and t 
are the root mean square value and the approx- 
imate correlation time of the specific stochastic 
force in one direction. 

The change in semi major axis a due to the 
effect of stochastic forces with diffusion coefficient 
D after time t is given by 



Ac 



= -VDt, 



(42) 



where uj is the mean mo tion of the moonlet 
( Rein and Papaloizou 20091) . Note that in this 
regime, the acceleration of the moonlet doesn't 
depend on its mass (as can be seen in Eq. (|4T|) ). 
Therefore, the migration rate and the diffusion co- 
efficient are independent of the mass of the moon- 
let ; this might be an observational indication for 
this migration regime. 

Assuming an initial semi-major axis of a = 
130 000 km and D - lO^^^m^s-^ as found in 
simulation SI, one can calculate the expected 
difference in semi major axis after one orbit 
due to stochastic forces which turns out to be 
Aa = 0.01m. For simulation S3, assuming 
D ^ 10^^^ m^ s^^, one finds Aa ~ 0.11m. For the 
simulation L2, assuming D ^ 10""'^^ m^ s~'^, one 
finds Aa = 10.5 m. These translate to random 



walks with standard deviation given as a func- 
tion of time by Aa — 0.27^t/(l year) m, Aa — 
2.7^t/{l year) m, and Aa = 270^t/(l year) m 
respectively. 

4.3. Discussion 

From the above, it can be seen that increasing 
the surface density by factors 4 — 5 changes the 
migration rate by two orders of magnitude. Thus, 
the results show clearly that the surface density S 
is much more important than in the regular, type 
I like migration, presented in Sect. [31 This can 
be easily understood with a toy model. The crit- 
i cal unstable w avelength A scales linearly with S 
(|Toomrelll96l . If we assume a fixed moonlet size, 
the ratio of moonlet size to A therefore changes 
with S. In the limit where A is much smaller than 
the moonlet radius, the stochastic forces are negli- 
gible as the density distribution is approximately 
homogeneous on the relevant scales. This is the 
case in simulation SI. In the other limit where A 
is larger than the moonlet, the moonlet undergoes 
a random walk that is similar to that of individual 
ring particles, as seen in simulation L2. 

The range in migration rates found in sim- 
ulations shows that over time scales of several 
years, the migration of a moonlet of mass m ~ 
10~-'^^Msaturn niay be dominated by a random 
walk in some situations (eg. those in simulations 
S3 and L2, the latter carried out with particles 
of radius 13 m). However, in regions of the rings 
where the surface density is small (eg. simula- 
tion SI), the moonlet may be in a regular, non- 
stochastic migration regime. In that case, the 
model from Sect. |3]can be applied. 

This dependence offers an exciting possibility 
to constrain the nature of the ring particles and 
the physical processes occurring in the rings by 
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Simulation 


Correlation time [s] 


Diffusion coefficient [m^ s '^J 

Dy 


Q 


L2 cutoff 
smooth 


5000 7000 
6000 9000 


9.61x10-^^ 14.97x10-^^ 
5.34x10-12 9.11x10-12 


4.1 


SI cutoff 
smooth 


2000 4000 
3000 8000 


1.92 X 10-" 2.26 X 10-" 
1.59x10-" 1.69x10-" 


7.2 


S3 cutoff 
smooth 


2000 6000 
4000 10000 


7.50 X 10-1'^ 30.24 x lO-^'' 
8.87x10-16 20.48 x 10-16 


2.5 



Table 2: Simulation results. The first column gives the name of the simulation, as defined in Table [T] The 
second and third columns give the correlation time in the x (radial) and y (aziniuthal) direction, respectively. 
The fourth and fifth columns list the diffusion coefficients. The sixth column is the Toomre Q parameter, as 
measured in the simulation. 



measuring the migration of moonlets. But note 
that because regular migration gives a decrease in 
the semi-major axis that is linear in time, provided 
it continues to operate, it will always ultimately 
dominate the behavior for large time because the 
spreading of the semi-major axis associated with 
stochastic migration increases only as the square 
root of time. 

At the present day, the number of observed mi- 
gration rates is not sufficient to draw a statisti- 
cally significant conclusion. However, the fact that 
the migration varies in rate and direction clearly 
favors the stochastic migration model presented 
in this section. Considering a migration rate of 
|Aa| = 100m in t = 1 year in Eq. P^ . one finds 
D ^ lA X 10-12 m^s-^. Taking |Aa| = 40m 
in t — 1 year gives I? = 2.2 x lO-i^m^s-'^. 
This is in the range obtained in the simulations, 
and tends to favor the case of simulation L2, and 
S 700 kgm-2 in the A ring. 

5. Conclusion 

In this paper we have calculated the differen- 
tial torque exerted on a moonlet by the outer and 
the inner disk with a smooth, flat surface den- 
sity profile. We performed both an accurate nu- 
merical integration and a second order analytical 
calculation. These approaches were found to be 
in excellent agreement where their domains of va- 
lidity overlap. The migration rate found in this 
case is proportional to the mass of the moonlet 
to the power 1/3. It is about —1 myr-i for a 
200 m radius moonlet in the A-ring. This is way 
too low to explain the observed migration of the 



propellers in Saturn's rings. Nonetheless, density 
fluctuations in the rings, due to their proximity 
to gravitational instability, can lead to stochastic 
torques on a moonlet, that may dominate on the 
timescale of the Cassini mission, to an extent that 
depends mainly on the surface density of the rings. 
These stochastic torques may account for the ob- 
servations. 

The possibility that the migration of propellers 
is induced by stochastic processes rather than by 
a regular type I like migration is therefore very 
exciting : this may help to infer the local surface 
density, and therefore the size of the ring parti- 
cles. Indeed, both quantities are linked through 
the optical depth, which is observationally well 
constrained. Our estimate of E ~ 700 kgm-^ is 
equivalent to '-^ 10 m size particles in the A ring 
for a single-sized population. This is in good 
agreement with available estimates from stellar oc- 
cu ltations. For the A ring. Voyager occultations 



(jZebker et al.lll985| ) flnd that the radius ra of par- 
ticles follows: 0.1 m (assumed) < Va < 11 m, 
with a power index ~ —3. For the 28 Sgr occul- 
tation ( French and NicholsonI [ioooh . a range Im 
< Ta < 20m is found, with a power index between 
—2.7 and —3, and an effective size (the single aver- 
age size accounting for t he fluctuations in photon 
count) being about 7m ( Cuzzi et al.ll2009[) . 

Fortunately the Cassini mission has been ex- 
tended to 2017. In the meantime, numerous ob- 
servations of the propellers will hopefully give a 
clear picture of their orbital evolution for a pe- 
riod of 10 years, representing about 10 000 orbits. 
This will allow us to test the hypothesis presented 
in this paper, and in particular to check whether 
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the propellers really are in stochastic migration, 
whereas first results seem to favor the stochastic 
hypothesis. 

Among the questions that still need to be 
addressed is why all the propellers seem to be 
gathered in the A ring, in places apparently de- 
void of density waves? Indeed the propellers 
seem gathered in a couple of narr ow radial bands 
of on ly about 1000 km width ([Tiscareno et al. 
20081 ). This is especially surprising since the 



the other hand, an estimate of the lifetime of a 
Pan size moon (~ 14 km in radius) against the 
today's cometary flux is provided by iDones et al 



A ring is densely populated by numerous density 
waves launched by the nearby small moons (At- 
las, Prometheus, Pandora, Janus, Epimetheus). 
Is there a systematic mechanism that would eject 
the propellers away from density waves ? Or does 
the present location of propellers just reflect the 
initial location of the parent body, assuming that 
the propeller population comprises the fragments 
resulting from the destruction of an ancient moon 
orbiting within the rings ? 

If the moonlets really undergo stochastic mi- 
gration, then Eq. may strongly constrain the 
age of the propellers, which can't be larger than 
the time needed to diffuse over Aa > 1000 km. 
Unfortunately, as long as D is unknown Eq. (H^ 
doesn't provide any useful information. However, 
considering that Aa is proportional to the square 
root of the time, and assuming that a moonlet 
migrates about 100 m in 1 year, one finds that 
Aa — 1000 km for t = 100 million years. Assum- 
ing |Aa| = 40 m in 1 year, we find that it requires 
~ 625 million years to diffuse over 1000 km. Note 
that for 200 m radius moonlets, the spreading due 
to stochastic migration equates to the contraction 
of the semi-major axes occurring as a result of 
smooth, regular migration after ~ lO'* yr, and 
then Aa « 10 km. Thus it would take about a 
million years to migrate through Aa — 1000 km 
in this case, largely through the action of the non- 
stochastic, regular migration process, if that can 
be assumed to operate smoothly and simultane- 
ously with the stochastic migration process. A 
100 m radius moonlet would migrate through only 
500 km in the same period, so that the smooth, 
regular migration process spreads a population of 
moonlets of various sizes over 1000 km in one to 
two million years. These could be the times since 
the catastrophic disruption of a small moon orbit- 
ing at 130 000 km from Saturn, that was broken 
into smaller moonlets by a meteoritic impact. On 



()2009l) and gives a range between 100 Myr and 
16 Gyr, depending on the size distribution of im- 
pactors. Therefore, the recent occurrence of such 
an event, about 4 to 4.5 billion years after solar 
system formation, is possible. Note also that an 
age of about 100 Myr is coherent with some esti- 
mates of Saturn's ring age despite of the lack of 
fully satisfactory exp lanation for their origin (see 
Charnoz et al.ll2009al for a review). 

We see that the question of propeller's migra- 
tion is inextricably linked to the issue of the ori- 
gin of Saturn's moons embedded in the rings, 
which is still a myst ery. iPorco et al.l (|2007l ) and 
Charnoz et al.l (j2007l ) have jointly proposed that 
small moons embedded in the rings could be ag- 
gregates of material on an initial shard denser than 
ice. When destroyed by meteoritic bombardment, 
these could release dense chunks of material that 
could explain the origin of the propellers. How- 
ever, the origin of Saturn's ring syst em is still 
a mat ter of debate ( HarrisI Il984: Charnoz et al. 
2009bl ). Knowledge of the age of the propellers 



could provide important constraints on the age of 
the main ring system and its embedded moons, as 
there are strong indications that these could have 
about the same age, provided t hese moonlet s hide 
a dense shard (Charnoz et al. 2007: Porco etal 



20071 ). Understanding the migration rate of the 
propellers is therefore an important piece of this 
puzzle. 
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and the organisers of the "Dynamics of Discs and 
Planets" workshop at the Isaac Newton Institute 
in Cambridge where these took place, as well as M. 
Tiscareno for providing us with migration rates. 
Hanno Rein was supported by an Isaac Newton 
Studentship, STFC, and St John's College, Cam- 
bridge. 
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A. Evaluation of the Green's function 

Here we evaluate the Green's function defined by Eq. (fT5|) as 



n=co I . „ ^ 



n— — oo 



To perform the summation we use the general result that if for a general periodic function 

n— oo 

K«)exp(m/3t), (A2) 

Tl — — oo 

with, period 2^/(3^ and h{n) being defined as an integrable function, we set 

F{t) = — / h{n) exp{inPT)dn, (A3) 



27r _ 
then 

n=oo 

.g(r) = 27r ^ F(T + 27m/^). (A4) 



We set 



Then the integral Eq. (jA3|) defining F{t) is readily performed by contour integration with the result that 
for t > 0, 

Fir) = -M-ir/2)sHu;,r) ^ ^^^^ 



otherwise F{t) = 0. Here uj^ — -^/TP^^t^/^. Using the above to evaluate the sum Eq. (|A4p as a geometric 
progression yields ^(t) = G(r) for < t < 2n/l3 as 

^ exp(-7r/2) sin(^-yr) - cxp(-77r//3) sin(^-y(T ~ 27r//3)) 

a;^^[l + exp(-277r//3)-2exp(-77r//3)cos(27rw^//3)] ' ^ ' 

the function is determined elsewhere by its periodicity with period 2ii/ (3. 
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